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Abstract 

We consider the if-body local correlations in the (repulsive) ID Bose gas for general 
K, both at finite size and in the thermodynamic limit. Concerning the latter we develop 
a multiple integral formula which applies for arbitrary states of the system with a smooth 
distribution of Bethe roots, including the ground state and finite temperature Gibbs-states. 
In the cases K < 4 we perform the explicit factorization of the multiple integral. In the case 
of K = 3 we obtain the recent result of Kormos et.al., whereas our formula for K — 4 is new. 
Numerical results are presented as well. 

1 Introduction 

The delta-function interacting ID Bose gas (also known as the Lieb-Liniger model or 
the Quantum Nonlinear Schrodinger equation) is one of the oldest and most important 
integrable models. Its study goes back to the papers [1, 2] where it was shown that the 
spectrum can be obtained by the Bethe Ansatz [3]. The thermodynamical properties of the 
model were determined in [4] using the method nowadays known as the Thermodynamical 
Bethe Ansatz (TBA). After these seminal papers tremendous effort was devoted to the 
calculation of correlation functions using various approaches [5, 6, 7, 8, 9, 10, 11]. One of the 
most important recent results is the exact determination of the long-distance behaviour of 
correlations [12, 13, 14, 15]. 

Apart from purely academic interest, the study of the ID Bose gas was spurred by the 
recent success of experiments with cold atoms in quasi one-dimensional traps [16, 17, 18, 
19, 20, 21, 22]. A remarkable result was presented in [19], where the authors managed to 
measure exact predictions of the TBA (for further developments and open questions see 
[23]). In experimental situations the local correlations are of special interest, for example the 
three-body local correlation is related to the rate of particle loss [24, 25, 22] and to the third 
moment of the density fluctuations [26, 27]. Moreover, even the four-body correlations might 
be accessible to experiment, as it was recently demonstrated in a 3D experiment [21]. 

Concerning the general if-body local correlations (for the precise definition see the main 
text) there has been considerable theoretical progress, too. The K = 1 case is simply given by 
the (linear) density of particles, whereas the K = 2 case was related to the thermodynamical 
quantities of the model in [28] . Concerning the higher-body cases small-coupling and large- 
coupling expansions were performed in [28, 29], whereas the exact ground state value of 
the three-body correlation was calculated in [30]. A new approach was initiated in [31, 32], 
where an infinite integral series (also called the LeClair-Mussardo or LM series) was derived 
using a special non-relativistic limit of the sinh-Gordon model. The LM series applies for 
any K and arbitrary temperature, including the ground state, and it can be considered as an 
effective large-coupling expansion of the quantity in question. The papers [33, 34] considered 
the relation between the LM series and previous form factor calculations with the Algebraic 



Bethe Ansatz (ABA) ; in [34] it was shown that the LM series can be understood and proven 
within the ABA. However, there was one crucial problem: there were no explicit and general 
results available for the form factors entering the LM series; the numerical results in [31, 32] 
were obtained using a truncation of the full series. 

The important task of the exact summation of the LM series was performed for the first 
time in the recent article [35], where the authors evaluated the three-body correlation based 
on a well-supported conjecture for the corresponding form factors. To our best knowledge 
this is the first time that an exact, explicit and compact result was given for a non-trivial 
correlation of the ID Bose gas, valid at arbitrary couplings and temperatures. 

In the present work we contribute to the calculation of the if-body correlators using a 
different approach. Our strategy is the following. First we consider a related physical quantity 
(the so-called "emptiness formation probability") on a generic XXZ spin chain and show that 
a special scaling limit of the spin chain [11, 36] yields the desired correlations in the Bose gas 
(Section 3). The matrix elements of the operator on the spin chain are calculated in Section 
4 borrowing results from the works [37, 38]. We then perform the scaling limit towards the 
Bose gas in Section 5, this way we obtain the form factors in a finite volume, with a finite 
number of particles. Finally, the thermodynamic limit is performed in the Bose gas (Section 
6) leading to the multiple integral (6.7), which is the main result of this work (see (6.10) for 
the dimensionless form). 

In principle the multiple integrals could be evaluated for any K, but in practice this 
becomes more and more difficult with increasing K, therefore it is desirable to derive more 
compact results. In Section 7 we show how the factorize the multiple integral in the cases 
K < 4. The results are the expressions (7.3), (7.10) and (7.12). In subsection 7.5 we also 
present examples of the numerical results. 

Finally in Section 8.3 we determine all form factors entering a modified form of the LM 
series, making it an explicit integral series for the i^-body local correlations. 



2 The Lieb Liniger model 

The second quantized form of the Hamiltonian is 

H LL = [ dx(d x ^d x ^> + c^^-W) . (2.1) 
J o 

Here L is the size of the system, periodic boundary conditions are understood and ^(x,t) 
and W{x,£) are canonical non-relativistic Bose fields satisfying 

[*{x,t),ti\v,t)]=8(x-y). (2.2) 

We used the conventions m = 1/2 and h = 1 and c > is the coupling constant. 

The eigenstates of the Hamiltonian (2.1) can be constructed using the Bethe Ansatz 
[1, 2, 7]. The TV-particle coordinate space wave function is given by 



7=E ex PvE^ (^ \ n 

V { 3 ) j>k 



(Pp)j - (Pp)k - ice(xj - x k ) 



\ v'/'ki = — > ,<■*]' {> > ) | | V " (Vph -(Vp) k (2 ' 3) 

where e(x) is the sign function. 

Periodic boundary conditions force the quasi-momenta to be solutions of the Bethe Ansatz 
equations 

e*» L TT P] ~ Pk ~ tC = 1. (2.4) 
The energy and momentum of the multi-particle state is given by 

E N = J2P 2 J P N=Y,P3- 

3 3 
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The norm of the wave function (2.3) is [39, 40] 



N LL = I \xn\ 2 = II {P \ Pk)2 t/ x detg" (2.5) 
J (P 3 -Pk) 2 



with 

/ N 



gf k L = s jtk [L + J2 vfo -pi))- ^fe - pi) ( 2 - 6 ) 
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We will be interested in the matrix elements of the operators 

O K = (^(0)) K (V(0)f . 
In coordinate space the matrix elements are given by the integrals 
AM 



<>n\Ok\xn) = 



Kl(N-K)\ 

L 

dxi... dx N -K (/)* N (0,...,0,Xi,..., x N - K )xn(0, ■ ■ ■ , 0, x\, . . . , x N - K )- 



(2.7) 



The expectation value of Ok describes the probability to have K particles at the same point. 
It is useful to introduce the dimensionless quantities 

{Ok) 

9K - —7^, 



where n = N/L is the particle density. It can be shown by scaling arguments that in the 
thermodynamic limit gx only depends on the dimensionless parameters 

c T 
1 = ~ r = — , 

n n z 

where T is the temperature (we used the convention ks = 1 for the Boltzmann constant). 

In principle the form factors (2.7) could be obtained by performing the integrals in coor- 
dinate space, but this becomes increasingly complicated with growing N. Note also that the 
Algebraic Bethe Ansatz for the Bose gas does not lead to simple results either: the action of 
the field operators on Bethe states can be evaluated easily, but afterwards one would have 
to compute scalar products of Bethe states with a reduced set oi N — K particles, neither of 
which are on-shell, and there is no good formula for the scalar products of such states. One 
way out of these problems is to consider a related quantity (the "emptiness formation prob- 
ability") on the XXZ spin chain, where there are methods available to compute its matrix 
elements. 

3 The XXZ chain and its special scaling limit 

The XXZ spin chain with M sites and periodic boundary conditions is given by the 
following Hamiltonian: 

M M 
H = JJ2( S * S ?+i + S J S j+i + A (S*S* +1 1/4)) + hJ2S?- (3.1) 

3=1 3=1 

This model is also solvable by the Bethe Ansatz [3, 41, 42, 43]. The A^-particle eigenstates 
are given by 

L L 

■m £ 4>N{\\yi,...,yN)a- x ...a- N \Q). (3.2) 



3/1=1 yi=N 
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Here |0) is the reference state with all spins up and ijj are the positions of the down spins. 
The amplitudes are 

, s 1 \- tt sinh((-PA) m - (VX) n + e(y n - y m )rj) 1 \ > 

^ (A|y) = 7M-^ U sinh((^A) m - (v\) n ) wnmuyti 

V V£a N l<m<n<N U m V w Z=l 



where 



(3.3) 

The parameter 77 is related to the anisotropy: 

A = cosh i]. 

In (3.4) we introduced inhomogeneities £j for the sites of the spin chain; they will be used as 
a technical tool to obtain the form factors in section 4. The physical limit consists of setting 
all £j — » 77/2. The expression (3.3) is a seemingly over-complicated way to write down the 
wave function, because it is valid at arbitrary values of the variables yj and not only in the 
region y\ < • • • < yu . We used this form to have an exact agreement with the conventions 
used in (2.3). 

The Bethc equations follow from the periodicity of the wave function and they read 

jj $M^±^ = !, (3.5) 
AA smh(A J — A fe — 77) 

where 

^) - n .y^v (3 - 6) 

aa smh(A - & + Tj) 
In the normalization (3.2)-(3.3) the norm of the wave function is given by 



^ MZ = E-E \<Kjn,...,vs)\ a = 

Vl Vn 

(-sinhr?)-^ J]/(A i) A fe )/(A fc ,A,) x det£ 

j<k 



xxz 



(3.7) 



with 



G * xz = ^ fc ( ^ + E - A;) I - ^™ (A ^ - Afc) - (3 - 8) 



The kernel ip xxz is given by 



xa-Z/,n -sinhT? 

smn(u + f)l Z) smn(u — 77/2) 

One-particle momenta and energies are given by the formulas 

^(X) = sinh(A + t?/2) = j sinh 2 /? _ ^ 

sinh(A — 77/2) cos(2A) — cosh?; 

3.1 Towards the Lieb-Liniger model 

There is a special scaling limit of the XXZ chain which yields the physical quantities of 
the Lieb-Liniger model [44, 36, 11]. In order to obtain the Bose gas in a finite volume L one 
has to set 

c 

77 = i~K — is M = —L 
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and let e — > (here c is the coupling constant of the Bose gas) . The number N of the magnons 
has to be kept fixed and the rapidities of the particles have to be scaled as 

After the limiting procedure the magnons can be identified as the particles of the Bose gas 
with rapidity pj . It can be shown that under an appropriate scaling of the parameters J and 
h 

f \ \ T sinh 2 T] 2 

e(A) = J 7— — n —> p — u, 

y ' cos(2A) - cosh 77 1 p ' 

where is the chemical potential in the Bose gas. However, this will be not needed in the 
following; we will consider the Bethe wave functions and the form factors of local operators. 
In the following we assume that the homogeneous limit £j — » 77/2 is performed first on the 
spin chain, and the limit towards the Bose gas is taken afterwards. 
Taking the scaling limit of the Bethe equations (3.5) results in 



For the sake of simplicity we only consider even chains so that no twist appears in the Bethe 
equations. 

2 

The limiting form of the Bethe wave function can be taken by setting Xj — ^rUj and 
keeping Xj finite, which will correspond to the position of the particles of the Bose gas. The 
wave function then reads 

* N (x\v) = ^ J] _ \[F((Pv) hXl ), (3.11) 

V Ve<J N m>n K > m v > n 1=1 

where 

F(p,x) = e- ivx {-l) y . (3.12) 

Apart from factors of (—1) the above expression is equal to the complex conjugate of the 
Bethe wave function (2.3). It can be argued that the factors of (— l) Vj don't affect the 
calculation of form factors of local operators. Indeed, for any coordinate space calculation 
one has to take the product of two wave functions with the down spins placed at prescribed 
positions. Depending on the operator in question an overall factor of (—1) may remain, but 
the position dependent factors of (—1)* always cancel. For the operators considered in this 
paper every such factor cancels, therefore they will be neglected in the following. 

Due to the relation between y and x it is expected that the norm of the wave function 
behaves as 

M xxz -> M LL . (3.13) 

Comparing the formulas (2.5) and (3.7) we obtain the same scaling: the Gaudin determinants 
behave as 

/ ■ \ N 

detg xxz -> (j) do-.; 



and the prefactors contribute an extra (is) 



-N 



3.2 The emptiness formation probability 

We are interested in the local operators E" acting on site j with matrix elements 



(e^~) =5 k , a 6i, p . 



H 

In particular we consider the composite operator 

s K =E^-E^~ ...E--. (3.14) 
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When sandwiched between two states, this operator forces K particles to occupy the first 
K sites. The expectation value of Sk (or sometimes its spin reverse) is called the "emptiness 
formation probability". 

We will show that the operator sk scales to Ok in the limiting procedure. In the coor- 
dinate Bethe Ansatz its TV-particle form factors are given by 

<{A}MM) 



K\{N-K)\ 

m (3.15) 

^2 0ArC*|l> 2 > ■ • • ' K > 2/1 > • • • ) VN-K^NifAh 2 > • • • ) K > Vl) • ■ • i UN-k)- 

yi,...,y N - K =K+l 

This formula has to be compared to (2.7). It is easy to see that the scaling limit of the Bethe 
wave functions works even if a fixed number of particles are placed on the first few sites: 

<f> N (fj,\l,2,...,K, y 1 ,...,y N - K ) -> Xn(p\0, 0, . . . , 0, x x , . . . , x N - K )*. 
Therefore the un-normalized form factor will behave as 

2 x M—K 

£ -) (MMM) {{ P }\o K \{k})\ 

where it is understood that 

c c 

-\. ' -/'; ' h \r 

For the normalized form factors we get 

' £ 2 v a '({a}i^im> {{ P }\o K mr 



c 



(3.16) 



Our strategy is to obtain explicit determinant formulas for the matrix elements (3.15) and 
to take the scaling limit according to (3.16). 



4 Form factors in the XXZ chain 



In this section we compute explicit determinant formulas for the matrix elements (3.15) 
in the framework of Algebraic Bethe Ansatz (ABA). Mostly we will use the results of the 
papers [37, 38]; the only difference between the present approach and the traditional meth- 
ods is that here the homogeneous limit £j — > rj/2 is taken explicitly before performing the 
thermodynamic limit or the limit towards the Bose gas. 

The central object in ABA is the monodromy matrix, a 2 x 2 matrix in the so-called 
auxiliary space with operator valued entries which act on the Hilbert space of the spin chain: 



T(u) 



A(u) 
C(u) 



B(u) 
D(u) 



It is built from the so-called local L-matrices: 



T(u) = L M (u) . ..Li(u), 

where 

Lj(u) =R € j(u-^ j ). 

Here j refers to the quantum space of the spin at site j and refers to the auxiliary space 
and the parameters £j are identical to the inhomogeneities already introduced in (3.4). The 
operator R(u) is the R-matrix of the XXZ type: 

/ sinh(it + rj) \ 



R{u) 



sinh(u + 77) 



sinh(u) sinh(?7) 
sinh(?7) sinh(u) 



(4.1) 



sinh(u + i]) J 
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The trace of the monodromy matrix is called the transfer matrix: 

t(u) = A(u)+D{u). 

In the homogeneous limit £j = 77/2 it is related to the Hamiltonian (3.1) at h = as 

d . . 

+ const. 



H ~ ±r(u) 
du 



--ri/2 

The normalization (4.1) of the i?-matrix results in the following vacuum eigenvalues: 

A(«)|0) = |0) D(u)|0)=d(u)|0), 

with d{\) given by (3.6). 

In the framework of ABA the Bethe states are 

N N 

(Oin^A,) and JJS(A*j)|0>. 

i=i 3=1 

They are eigenstates of the transfer matrix if the rapidities satisfy the Bethe equations (3.5). 
Apart from an overall normalization factor they are identical to the states given by the 
coordinate wave functions (3.2). 

In order to compute the matrix elements of sk hi the framework of ABA the local 
operators have to be expressed in terms of the entries of the transfer matrix. This problem 
was solved in [38, 45] leading to the following theorem: 

j-i m 
Ef = JJ (A + !>)(&) xT""&)x JJ (A + £>)(&)• (4-2) 

fc=l k=j+l 



Applying this formula to operators E- on neighbouring sites one gets 

L 

8 k = D(t 1 )Dfa)...Dfa) JJ (A + D)&). (4.3) 

l=k+l 

Evaluated on Bethe states equation (4.3) yields 

N N N N k 

(0| JJ C(Xj) s k JJ B(^)\o) = (0| JJ C(\j) £>(a)£>(6) • • ■ IJ B(^)\0) x JJ 
j=i j=i 3=1 j=i 7 -=i 

(4.4) 

with i(it,{//}) being the corresponding eigenvalue of the transfer matrix. Evaluated at the 
inhomogeneities it reads 

N 



fc=i 

In (4.4) we also used the fact that 



JJ(a + j d)(6) = i- 



z=i 



The action of multiple D operators on the dual state results in [37] 



N 



(0\l[C{\ j )D(Z 1 )D&)...D(l; k ) = 

i=i 

2^ ff f \ ■ r7T+ r+T x dett A^,0) x JJ/(A+,AJ: 

{a + ma-} smh (& ~ G)smh(AJ - A+) ^ 

|{A+}|=fc 

JJrf(A^)<0|JJC(Ao)C(6).-.C(4) 
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with 



f(A,0 = 



sinh?7 



sinh(A - £) sinh(A - £ + 77) 
The scalar product of an arbitrary state and a Bethe state is [46] 

Uj,k sinh(^- - \k+v) 



inc(A,)n^-)io)= 



where 



Yl j<k sinh(Aj - Afe) sinh(/_t fc - 



JV 



x det S, (4.6) 



S'jfe = t(fXj,X k ) - d(X k )t(Xk, Hj) Yl 



sinh(A fc - fu + 77) 
1 sinh(A fe - /i; - 77) ' 

Specializing this to the present case 

<oi n c (k MM ■ . • c(6) n = 
° j 

UjM sjnhQxj - A~ + 77) ]1 3 Q sinh(/j J - £ Q + 77) 

ELk* sinh(/J fc - Mj) rij</c sinh(AT - A ,7) ]\ J<k sinh(^- - £ fe ) ^ fc sinhfe- - A fe ) 
x det U. 



Here 



Un=t{jii,Zi) for l = l...K 



N 



„ = i smh(A z - /J - 77) 

Therefore, the form factor of the inhomogeneous chain is given by 



otherwise. 



(4.7) 



<{A}klM> 



n 



1 



Uo, P smh(A+ - Co + t?) 



f=\ X {A+ ^A- } Hi<i Binh(6 - 0) sinh(At - A+) ^ 



Il/(A^A-)x J]d(A+) 



|{A+}| = fc 

lli.fc sinh^ - X k + 77) n 3 > sinh(/ij - £ G + 77) 



(4.8) 



H j<k sinh(/j fe - /jj) n j<fe sinh(A j - X k ) ]J j<k sinh(£j - £ fc ) ]lj,fe sinhfo - A fc ) 
x dett(A^,£j) x detU. 

We now perform the homogeneous limit £j — > 77/2 following the method of [47]. For the 
matrix A/ we get 

det M 1 



lim 



n^sinh^-a) nri 1 



■ det H 



with 



H 



.1 1 



9£ 



1-1 



t(x k ,o 



,1 



=r,/2- 



It is advantageous to use the form 



*(A,0 = 



cosh(A — £) cosh(A — £ + 77) 
sinh(A - £) sinh(A - £ + 77) ' 

Taking the derivatives one is free to replace [48] 



l-X 



t(x,0 



(-iy-\i-iy. 



( cosh(A — £) \ ' / cosh(A — £ + 77) 



\sinh(A — £)/ \ sinh(A — £ + 77) 



S 



The same steps can be performed for the corresponding elements of the matrix U. Finally 
the homogeneous limit reads 



{A+}U{A-} lW Sllll H A j " A i J o,p i 
|{A+}|=K 

IIj,fc sinh(/ij - Afc + t?) n 3 sinh g (/Xj - rj/2) 

Uj<k smh(pj ~ Mfe) ]lj<fc sinh(AT - A") J]/ sinh K (? ? /2 - A;) 



(4.9) 



x det O det V. 



with 



and 



On 



( cosh(A+ - 77/2) \ / cosh(A+ + ? ? /2) ' 
^sinh(At-r7/2) ) ~ \ sinh(A+ + V /2) / 



(4.10) 



cosh(/ij — 77/2) \ ' / cosh(/x,- + 77/2) \ ' 
sinh(/ij - 77/2) / \ sinh(/ij + 77/2) / 



for l = l...K 



N 



(4.11) 



V*,K+i = t(fi j ,X l ) ~ d{\ )t{\ >Mj) n 



otherwise. 



sinh(A ; - Ho + i]) 
sinh(A ; ~ -n -ri) 

5 The scaling limit of the form factors 

Here we take the scaling limit of the formula (4.9) to obtain the matrix elements of Ok 
in the Bose gas. We substitute 

f] = VK — is A = — p [1 = -k. 

c ' c 

It is straightforward to calculate the limiting values of the prefactors, but the determinants 
need special care. The elements of O read 



The leading terms will be 



' sinh(A+ + is/2) \ 1 ( sinh(A+ - ie/2) \ ' 
cosh(A+ +ie/2) ) ~ \ y cosh(A+ - ie/2) J 



(5.1) 



O, 



which yields 

detO ->■ if!(«e) 



/r \ (if-l)if/2 

)*(£) det[(p+y-i]=K!( i£ : 



K 



(k— l)fc/2 



One can use the same expansion for the first if columns of the matrix V. 

To obtain the proper normalization note that in ABA the norm of the Bethe state scales 

as 



(oinc(A,-)n 5 ( A i)i°) 



N Ui<i((pj -pi) 



x det Q LL . 



Uj<l(P3-Pl) 2 

This differs from (2.5) by the overall factor of c N . 

Collecting all the factors and performing a complex conjugation we find 



{{ P }\o K \{k}) = c K - N {K\f Uf(p^Po)U e ~ iLp ' 

{ P +}u{p-} 0,1 1 

\{p+}\=K 

Ui,l( k 3 - Pi + ic ) 



(5.2) 



IIj<i(*i - k i)Ui <k (Pi -Pk) 



x det Z, 
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with 



Zji = (kj) 



for I = 1 . . . K 



N 



z jtK+l = tfapT) - e- iL vrt(p-, kj ) n^- fco + lc) 



0=1 (Pi - k ° - ic ) 



otherwise. 



(5.3) 



Here we used 



t(u) 



u(u + ic) 

Equation (5.2) refers to the normalization where the norms of the states \{p}} and \{k}) are 
given by (2.5). 

The result (5.2) is valid whenever the set {k} satisfies the Bethe equations. In the case 
when {p} is also a solution but different from {k} we obtain the form factors 



f% (m, {k}) = c k - n (k\) 2 n f(Po,pt) n x 

{ P +}u{p-} o,l 1 
|{p+}|=K 

Uj,l( k 3 'Pi + i0 ) 



(5.4) 



U^-kOU^iP 



x det V, 



j<kWj 



Pk. 



with 



for 1 = 1 ... K 

N 



Vj, K +i = t(k j ,p l ) + t{ Pl , kj )Y[ 



(Pi ~k + zc)(p ; -p - ic) 
aJ[ (Pi ~k - ic)[pl -p + ic) 



otherwise. 



(5.5) 



This is a new result of the present work. In the case of K = 1 eq. (5.4) yields an alternative 
representation for the form factors of the density operator, which were previously determined 
in [46, 8, 49]. 

To obtain the mean value of Ok we take the limit {p} —> {k} in (5.2) and divide by the 
norm (2.5) resulting in 



(<D K ) N = [K\f 

{p+}u{ P -} 

\{p+}\ = K 



n 



pf -pi 



j>l 



(pJ-pT) 2 + c 2 ) 



det U 
dot Q LL ' 



(5.6) 



The elements of T-L are given by 



(Pj) 1 - 1 for l = l...K 
gtf for l = K+l. 



.N. 



(5.7) 



Here it is understood that in both Q LL and H the ordering of the rapidities is given by 
{p} = {{p + }, {p~}}. The matrix T~L differs from Q LL only in those columns which belong to 
the subset 

In the case of K = 1 the above formula results in 

as it should. To prove this note that the sums of the columns of Q LL are equal to L in every 
row, therefore every T-L gives 

n = \g LL . 
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6 Expectation values in the thermodynamic limit 

In this section we evaluate the thermodynamic limit of (5.6). We consider a Bethe state 
|f2) in a large volume L with a large number of particles such that the particle density 
n = N/L is fixed. In the thermodynamic limit one defines the density of roots p^ r \p) and 
holes p( h \p) and the total density p(p) — p^ r \p) + p^ h \p). This latter function satisfies the 
Lieb-equation 

p(p) = ?T+ ( dq^(p- q)f(q)p(q). (6.1) 

Here 

p [r) (v) 

m - l$ (6.2) 

is a distribution function characterizing the state in question. In thermal equilibrium f(p) = 
(1 + e 6 ^) -1 where e(p) is the so-called pseudo-energy, which is a solution of the TBA 
equation 

e{p) = |j^(p ~ P') log(l + e-*'>). (6.3) 

At zero temperature we recover the ground-state distribution 

'«-{; S>a 

with A being the Fermi-rapidity. The particle number is always given by the formula 

N f 
n =-[= f(p)P(p)- 

We proceed to calculate the thermodynamic limit of (5.6) using the techniques of [50]. 
The ratio of determinants is calculated as 



detn 
detg LL 



The resulting matrix on the r.h.s. will be equal to the identity matrix except for those columns 
belonging to the set These elements can be evaluated by transforming the action of 

Q LL into an integral equation. This results in 

det H -i-r 1 , T 

-; — tttt = x- x det /, 

det G LL 11 2ttL P ( P +) 

where 

Here (u) is the solution of the linear integral equation 

h {l) (p) =p' + J°° ^ <p(p - <i)f(q)h {l) (<?)• (6.5) 

Note that in this normalization h^(p) = 2irp(p) and f{p)h^\p) = 2np( r \p). 

As a final step one integrates over the rapidities p~^ and in the thermodynamic limit one 
gets 

WW - k, J % . . . £ n , M n 57 i^ x m 

The prefactors are completely anti-symmetric therefore one can expand the determinant and 
write 

mew =<*!)»/ f ...^n s^ njw^ <«) 
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This formula is the main result of our paper. Its explicit factorization is performed in the 
cases K = 2, 3, 4 in section 7. 

For practical purposes it is useful to derive the dimensionless multiple integral for the 
quantity 

a - { ° k) 

which only depends on the dimensionless coupling constant 7 = c/n and the dimensionless 
version of the distribution functions f(p). We define 

9=- fiv) = f(p = ic)- 

In the finite temperature case f(q) = (1 + e^ 9 -*) -1 where e(q) is the solution of the dimen- 
sionless equation 

%) = T d J-- 2 —— log(l + e"^)), (6.8) 

r }-oo 2tt (q - q') 2 + 1) 

with 

/x T 
a = — t = —. 

T n 2 

Defining the dimensionless functions hS l > (q) as 

we find hS 1 ' (q) = c L h^ (p = qc). Thus the dimensionless multiple integral formula is expressed 
as 

« - {K ^ K It- ■ t n fe^FTi n /«,>*"-"<*>• («») 



6.1 The c 0+ limit 



We take the small-coupling limit of the dimensionful formula (6.6) by sending c — >■ and 
keeping n fixed. The limiting form of the kernel ip (u) is given by 

<p(u) -> 2nS(u). (6.11) 

Therefore the solution of the integral equations (6.5) is 



l-/(p) 

The determinant in (6.6) has the limiting value 

det/ ^ n iz7(p) 11^ -p») 



resulting in 



Note that 



therefore 



2=1 lj 



" '//'; /(Pj) 



(Oi) - n : 



27T 1 - /(ft) 

dp f(p) 

27Tl-/(p)' 
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as it should be for free bosons. 

The above calculation only applies if f(p) < 1, therefore the result is not valid for the 
ground state. In fact 

lim lim qk = 1. 

We have checked this property numerically for K < 4 (see subsection 7.5). To prove it 
analytically from (6.6) one has to carefully analyze the sub-leading corrections to the integral 
equation (6.5) with the weight function (6.4). We leave this problem for further research. 

6.2 The c ->■ oo limit 

Here we derive the leading term in the large coupling expansion of gx ■ In the c — > oo 
limit the kernel tp{p) is of order 1/c, therefore to leading order 



hW(p)=p 



and 

j>l 



Evaluating this formula for the ground state gives 

K(K-l) 



gK= 2K\^jJ X J 1 ■■■ dx K[[(Xj -Xi) 



3>l 

This result was already obtained in the papers [28, 29, 31, 32]. 



7 Factorization of the multiple integrals 

In this section we perform the factorization of the multiple integral formula (6.7) in the 
cases K = 2, 3, 4 1 . In order to keep the formulas as short as possible we will use the following 
notation: 

/#■•-/*/<,>•■■ 

Moreover we will suppress the dependence of the mean value on the state |f2) and we write 
simply {Ok)- The dependence on \Q) is carried by the functions f(p) and h^{jp). 

We found it more convenient to work with the dimensionful formula (6.7), because this 
way non-trivial checks of dimensional analysis can be performed at each step of the calcula- 
tion. The dimensionless formulas can be obtained as explained in 7.5. 

The main idea behind the factorization procedure is simple: at each step the number of 
the integrals can be reduced by one using the integral equation (6.5), whenever the prefactors 
are such that the corresponding variable is present only in one denominator 

7 \ 2 , 2 = 7r<p(pj ~ Pk)- 

Except from the case K = 2 this is not the case, instead the prefactors have to be divided into 
several terms, in each of which one of the integrals can be performed. This is a non-trivial 
task with growing complexity as K increases. In the following we present a case-by-case 
study up until K = 4. 

x The case K = 1 is trivial and it simply yields the particle density as it should. 
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7.1 K = 2 

One has 



(0 2 ) = 4 



dpidp 2 



P2 



(Pi -P2) 2 



Pi 



' (Pi - P2) 2 + c- 



■h^( pl )h^\p 2 



In the first term one integrates over p\ first, in the second term over p 2 first. Using the 
integral equation (6.5) one gets 



(O2) = ~ 



dp 2 p 2 /i (1) (P2)(^ (0) (P2) - 1) - / dpi Pi^ (0) (pi)(/i (1) (Pi) 



■Pi 



- / dp(j> 2 hw(p)-phw(p)). 



(7.1) 



This is in agreement with formula (10) of [35]. As it was already explained in [31, 32, 35], 
in the finite temperature case (7.1) agrees with the result obtained from the Hellmann- 
Feynmann theorem. We also note that (7.1) can be proven for arbitrary weight function 
f(p) using the Hellmann-Feynmann for a single state; one has to repeat the arguments of 
Appendix D of [34]. 

For future use we define 



{n,m} = J dpp n h^( P ) = J 0/(p) p n h^(p) 

It can be shown using the iterative solution to (6.5) that in general 

{n, m} — {m, n}. 

Using this notation we write 

2 



(7.2) 



(7.3) 



<(9 2 > = -({0,2}-{l,l})). 
Also, it is useful to derive a general formula which will be used often: 

dxdy X "~f +c2 {h^ {x)h,™ (y) - hW (y)h^ (x)) = 1 ({7, a + 0} - {f3, a + 7 }) . (7.4) 



7.2 K = 3 

The mean value is given by 



(0 3 ) =36 f #!...dp 3 n 

J 3>l 



Pj - Pi 



(Pi ~Pi) 2 + £■ 



■h^{ Pl )h^{p 2 )h^(p 3 ). 



(7.5) 



In this form neither of the integrals can be done directly. Instead, one has to divide the 
prefactors into several terms such that in each of them one integral can be performed. One 
way to do this is as follows. We define 



D(x,y,z) 



y-x 



(x — y) 2 + c 2 (y — z) 2 + c 2 



Then we find the identity 



n 



pj - pi 



3>l 
1 



(Pj ~Pi) 2 + c 2 



(7.6) 



[D(x, y, z) + D{y, z, x) + D(z, x, y) - D{y, x, z) - D(x, z, y) - D(z, y, x)} 



For simplicity we used the variables x, y, z on the r.h.s. instead of pi,p 2 ,P3- Equation (7.6) 
can be proven as follows. The r.h.s. is a completely anti-symmetric function of the variables 



14 



(7.7) 



x, y, z and it has exactly the same poles as the function on the l.h.s., therefore it has to be 
equal to the l.h.s. multiplied by a symmetric polynomial. By power counting it is shown that 
this polynomial is a pure number. This number is found to be 1 by simple manipulations 
after sending c — > on both sides. 

Substituting (7.6) into (7.5), performing one integral in each term using (6.5), changing 
variables accordingly and observing the cancellation of the terms including three h functions 
we find 

(0 3 ) = - [ dxdy - y ~ X x W(hi(x)ha{y) - h {x)h 1 {y))+ 
c J \x — yy + c z 

y(h (x)h 2 (y) - h 2 (x)h Q (y)) + 

{hiWtniy) - hi(x)h2{y))]. 

Using (7.4) the last line of (7.7) yields 

- / dM V ( V ~u± 2 {h2{x)h l {y)-h l {x)h 2 {y)) = ^({2,2}-{\^} 
c J [x — yY + c z c z V 

For the second line of (7.7) we can drop the term proportional to xy to find 

~ / did y {x _ y)2 + C 2 (ho(x)h 2 (y) - h 2 (x)h (y)) = 1 ({0, 4} - {2, 2} 
The first line of (7.7) is more involved. We have to compute 

>! f dxdy {y ~ X){ f 2 ± p (hi(aOfto(y) ~ M*)Mv))- (7-8) 
[x — yy + cr 



c 

We write 



(y — x)(y 2 + x 2 ) y~x c 2 x — y 2 y 3 — x 3 



(7.9) 



(x — y) 2 + c 2 3 3 (x — y) 2 + c 2 3 (x — y) 2 + c 2 

The first term in (7.9) gives 

^({0,1} 2 -{0,0}{1,1}). 

The second and third terms in (7.9) can be evaluated using (7.4). 
Putting everything together 

(0 3 ) =1 ( - 4{1, 3} + 3{2, 2} + {0, 4}) + ({0, 2} - {1, 1}) + ~ ({0, l} 2 - {0, 0}{1, 1}) . 

(7.10) 

This is in accordance with formula (11) of [35]. 

7.3 K = 4 

We define 

1 

Di{x, y, z, u) 



{{x - y) 2 + c 2 )((y - z) 2 + c 2 )((z - u) 2 + c 2 )' 
Then we find 

This equation can be proven through the same steps as in the case of (7.6). 
One has to evaluate 

(0 4 H 24 2 / djk...dpt]l ff Pi ~* 2 , h^{ Pl )h^{p 2 )h^{p 3 )h^{p i ) = 
J (CPj - Pi) + c ) 

r i 4 

= 48 / dp x . ..#477 no ^77 si ^T77 ^ 57 V (-1) [P] TT ^ 

J Pi-P2 2 +c 2 )((p 2 - p 3 ) 2 + c 2 )((P3 - Pi) 2 + c 2 \\ 
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Performing the integrals over pi and P4 and using the symmetries one gets 



There are in total 6 different combinations and we treat them one by one. We define 

12 



o=l 



with i^o representing one of the six combinations. The first three cases are evaluated easily: 
" ' : (^ 2 ) (x)^ 3 ) (y) - h& {x)h V) {y) ) = I ( {2 , 4} - {3, 3}) 



K\ = I dxdy 



(x — y) 2 + c 2 



K2 = I dM y {x y . y) 2 X +C 2 ( ft(3) fa) - hil) (^ (3) (y)) = ^ ({3, 3} - {1, 5}) 

(h™ (x)h^ (y) - {x)hW (y)) = \ ({1, 5} - {2, 4}) 



,,3 ^3 



K* = 



(a; — y) 2 + c 2 K c ' 

The remaining three case are more complicated because we have to separate factors of the 
form 



x a yP — x^y a 
(x — y) 2 + c 2 



a, (3 > 0. 



The next case is 



Here we use 



xy 2 — x 2 y x — y c 2 y — x 



1 y 3, - x 3 



(x — y) 2 + c 2 3 3 (x — y) 2 + c 2 3 (x — y) 2 + c 2 

leading to 

Ka = \ ({0, 1}{0, 3} - {0, 0}{1, 3}) + I ({0, 4} - {1, 3}) + 1 ({0, 6} - {3, 3}) . 



The next case is 



K, 



= J dM y { ^ y) 2 + c 2 (h {2 \x)h^\y) h«>\x)hW(y)). 



Here we use 



xy A — x y 1 



(x - y) 2 +c 2 2 



9 2 A A 

2 2 2 x - y y - x 

x-y- c z -, NO , n + 



(x — y) 2 + c 2 (x — y) 2 + c 2 
which gives 

K 5 = ({2, 2}{0, 0} - {0, 2} 2 ) + I ({2, 2} - {0, 4}) + 1 ({2, 4} - {0, 6}) . 
Finally, the last term is 



„2„,3 „,2_3 



K, 



= J dM \ x Zy/+c 2 ^ (0) {x)h(1) {V) - '' ! 1 ' ' ' ! ! : ! " ! ' 
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Here we write 



x 2 y 3 — x 3 y 2 1 
(x — y) 2 + c 2 5 
2c 4 



2c 2 

— {x - y) + (x 3 - y 3 ) + 2{x 2 y - y 2 x)+ 



y-x 



5c 2 



y 3 — x 3 



y 5 - x 5 



3 (x — y) 2 + c 2 3 (x — y) 2 + c 2 [x — y) 2 + c 2 



leading to 



K 6 — ({0, l} 2 - {0, 0}{1, 1}) + - ({0, 3}{0, 1} - {1, 3}{0, 0} 

~ ({0, 2}{1, 1} - {0, 1}{1, 2}) + ^ ({0, 2} - {1, 1}) + I ({0, 4} - {1, 3} 



(7.12) 



1 ({0,6} -{1,5} 
Putting everything together 

(Oi) = [8c 3 ({0, l} 2 - {0, 0}{1, 1}) + 32c({0, 1}{0, 3} - {0, 0}{1, 3} 
24c({0, 2}{1, 1} - {0, 1}{1, 2}) + 30c({0, 0}{2, 2} - {0, 2} 2 
4c 4 ({0, 2} - {1, 1}) + 5c 2 ({0, 4} - 4{1, 3} + 3{2, 2} 
{0, 6} - 6{1, 5} + 15{2, 4} - 10{3, 3} 
This a new result of the present work. 

7.4 Galilei invariance 

The expectation values (Ok) are Galilei invariant, and it is useful to check this property 
in our final formulas. This constitutes a highly non-trivial check, as it was already remarked 
in [35]. 

In our calculations we did not restrict ourselves to symmetric distributions, the weight 
functions f(p) can be arbitrary. Therefore, to check Galilei invariance it is enough to consider 
an infinitesimal boost b. This boost yields the following infinitesimal transformations: 



p — > p + b 



h u \p) 



h {:i) (p) + bjh U - 1) (p). 



It is then readily seen that (6.7) is invariant due to the anti-symmetry of the prefactors. 

We also performed the check on our factorizcd formulas. The transformation rules for the 
quantities {a, /?} are 



{a, (3} 



{a, (3} + b a{a -1,(3} + f3{a, ft - 1} 



Using this rule we have checked that the variation of equations (7.3), (7.10) and (7.12) indeed 
vanishes. 

7.5 Dimensionless formulas and numerical results 

The dimensionless versions of formulas (7.3), (7.10) and (7.12) are obtained simply by 
setting c = 1, multiplying with an overall factor of 7 , and replacing 

{n,m} -> Jdqq n U m \q), (7.13) 

To numerically evaluate the factorized formulas the following steps have to be performed: 
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logio 7 



Figure 1: The ground state values of the iT-body local correlations for K < 4 as a function of the 
dimensionless coupling j (g% = 1 by definition). The exact values are represented by the solid 
lines, whereas the dashed lines show the empirical formula (7.14). 



• Solve the TBA equation (6.8) iteratively. The parameter a can be fixed by requiring 

51 = 7 {0,0} = 1 

• Solve the linear integral equations (6.9) for h^(q). 

• Evaluate (7.3), (7.10) and (7.12). 

We peformed this procedure for a wide range of the parameters 7 and r. The quantity (74 
shows the same qualitative behaviour as 52 and (73 [35]: it is an increasing function of r and 
a decreasing function of 7, with the limiting values given by 

lim 34 = lim 34 = 4! = 24 lim lim (74 = 1. 

7— >0 t— >oo 7— >0 t- >0 

To demonstrate the numerical results we present the ground state values of £72, ff3 and g 4 
in Fig. 1, whereas the temperature dependence of (74 is shown in Fig. 2 for the intermediate 
couplings 7 = 0.1, 7 = 1 and 7 = 10. 

At T = the first term in the small coupling expansion of gx is given by [29] : 

We found that the empirical formula 

9KaaL j.m^i^) (7 , 4) 



gives a surprisingly good approximation and can be used for practical purposes even at 7 ~ 1. 
The predictions of (7.14) are also plotted in Fig. 1. It is expected that (7.14) holds with a 
good approximation even for higher K. 

It would be useful to compare the exact numerical values for g% and 54 case to the various 
approximations available in the literature [28, 29] including the large-coupling expansion both 
at zero and finite temperatures. This is out of the scope of the present work and is left for 
further research. 
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-1 1 2 3 4 5 

lo S io r 

Figure 2: The quantity as a function of the dimensionless temperature r for intermediate 
couplings. In the r — > limit the three curves approach small (but non-zero) values which are 
shown in Figure 1. 

8 Mean values in the LeClair-Mussardo formalism 

In this section we elaborate on the LeClair-Mussardo formalism, which is an alternative 
approach to obtain expectation values of local operators leading to an infinite integral series 
[51, 52, 31, 32, 34]. For our present purposes the following form of the series is the most 
convenient [34]: 



where 




and f(p) is defined in (6.2). The form factors appearing in the above series are defined as 

f£> 1: . . . , PN ) = n , {p r\f+ 2 x + ^ {Pi», (8.2) 

where the form factor on the r.h.s. is given by (5.4). This prescription is also called the 
"symmetric evaluation of the diagonal limit"; note that this limit is different from the way 
we obtained the mean value (5.6) because in (5.4) the Bethe equations were substituted into 
the matrix element before taking the diagonal limit. Therefore the object in (8.2) does not 
depend on the volume L. Note also that the l.h.s. refers to a normalization where the norm 
of the Bethe state is given simply by the Gaudin-determinant Q LL . 
Alternatively (8.1) can be expressed as [34, 31, 32] 




Here c (pi, ■ ■ ■ ,Pn) are the so-called connected evaluations of the diagonal form factors. 
Their precise definition and the relation to s (pi, ■ ■ ■ ,Pn) can be found in [52, 34]. The 
series (8.3) was originally developed in [51] in the framework of integrable Quantum Field 
Theories. Later it was used in [31, 32, 35] to compute the quantities gx up to K = 3. 
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However, for higher K the results (8.1)-(8.3) are only formal because the form factors them- 
selves were not calculated previously. In [31, 32] a prescription was given of how to obtain 
the connected evaluation using a special non-relativistic limit of certain form factors of the 
sinh-Gordon model. However, the actual calculation becomes more and more demanding 
with higher K and N . 

We fill this gap here by calculating the explicit results for Fjy s (pi, . . . ,Pn) for arbitrary 
K and N: we take the symmetric diagonal limit of the form factor (5.4). Note that every 
singularity of the form factor is included in the matrix Z, and even the elements of Z can 
be evaluating easily. Taking the limit and multiplying with the prefactors we obtain 



F^ s ( Pl ,...,p m ) = (K[) 2 £ 

{p+}u{p-} 



n 

3>l 



Pj 'Pi 

{pt-ptf + c 2 ) 



x det Y. (8.4) 



\{p+}\=K 

The elements of Y are given by 

= (Pi) 1 " 1 if pi£{p + } 

(8.5) 



Note that 



Y = H 



L=0 



where % is the matrix defined in (5.7). 

With these results the series (8.1) can be considered an explicit representation of the 
mean value. 

It would be desirable to have a general recipe for the re-summation of the series, which 
would be an alternative way to obtain factorized formulas like (7.10) and (7.12). However, 
this is far from being easy. The simpler cases K — 1 and K = 2 were already calculated 
in [31, 32]. The highly non-trivial case of K — 3 was considered in [35], where the authors 
evaluated the series (8.3) (and obtained the result (7.10) for the first time) based on the 
following conjecture for the quantities Ffj 2 : 

F$t, c = -^^2^12^23 ■■ .(Pn-i,nPin(Pin ~Pi2 ~ p\z Pn~i,n)- (8-6) 

Here we check this formula in the first two cases. In the simplest case of N = 3 our formula 
(8.4) gives 

P 2 i 

F ls = Fi <c = 36 J[ 2 3 2 - 
j>i P i l + C 

This was already calculated in [31, 32] and is in agreement with (8.6). In the case of N = 4 
one has to use the following relation between the symmetric and connected evaluations [52]: 

FI s {pi,P2,P3,Pa) = Fl c (pi,p 2 ,P3,p4) +X]- F 3,c (Pj) X I ^Vjk 

J \k^j 

Here pj means that pj is not present among the arguments of the form factor. We used the 
program Mathematica to express F% c using the above relation and we found agreement with 
(8.6). This is a highly non-trivial check of the conjecture (8.6); a proof for arbitrary N is not 
known. 

Finally we note that in the simpler cases of K = 1 and K = 2 we evaluated (8.4) and 
found exact agreement with the corresponding formulas of Appendix D in [34] . 



2 In order to ensure compatibility with our normalizations we inserted a factor of 1/c 2 
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9 Conclusions 



We developed multiple integral formulas for the local correlations in the ID Bose gas. 
The final results for the expectation value (Ok) is given by equation (6.7), whereas the 
dimensionless formula for gx is given by (6.10). 

In section 7 we performed the explicit factorization of the multiple integrals in the cases 
K = 2,3,4; for K — 3 we obtained the recent result of [35] whereas our formula for K = 4 
is new. Our method of factorization relies only on the integral equation (6.5) defining the 
auxiliary functions entering the multiple integral. Therefore the process works for arbitrary 
distribution of Bethe roots and not only for the ground state or the finite temperature Gibbs 
states. 

The general recipe of how to perform the factorization for K > 4 is not known. The 
strategy is clear: at each step the prefactors have to be manipulated in such a way that the 
number of integrals can be reduced by one using the integral equation (6.5). We believe that 
this can always be done and it would be interesting to develop a general algorithm for this 
process. 

An alternative way to obtain the mean values (Ok) would be to take the thermodynamic 
limit on the XXZ spin chain first and to perform the scaling limit towards the Bose gas 
afterwards. The advantage of this approach would be that on the spin chain the factorization 
of the multiple integral formulas for the elements of the reduced density matrix is by now 
well-understood (see [53] and references therein). In fact we attempted to take the scaling 
limit of the factorized results of [54] concerning the emptiness formation probability. However, 
this turned out to be cumbersome already in the case K = 2. Thus it seems that the direct 
approach of the present paper is more advantageous, at least for the small values of K 
considered here. 

An other alternative way would be to sum up the LeClair-Mussardo series (8.1) or (8.3). 
The diagonal form factors entering (8.1) are given explicitly by (8.4), therefore the remaining 
task is purely combinatorial: one has to expand the sums of determinants appearing in (8.4) 
and put the resulting expression in a form which is amenable for re-summation. Again, this 
is a formidable problem, the solution of which is not yet known. 
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